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Abstract 

The Peierls stress for a [lll]-screw disfocatfon in bcc Tantalum is calculated using an embed- 
ded atom potential. More importantly, a method is presented which allows accurate calculations 
of the Peierls stress in the smallest periodic cells. This method can be easily applied to ab ini- 
tio calculations, where only the smallest unit cells capable of containing a dislocation can be 
convicntly used. The calculation specifically focuses on the case where the maximum resolved 
shear stress is along a {110}-plane. 
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1 Introduction 

Long, low-mobility [lll]-screw dislocations control low temperature plastic behavior in bcc metals. 
Unlike their fee counter parts, bcc metals violate Schmid behavior and have many active slip planes. 
The microscopic origins of such behavior is crucial in the understanding of their plastic behavior. 
Therefore detailed and accurate first principle calculations are invaluable. 

The Peierls stress, the zero temperature limit of the critical stress for slip, has been calculated 
in a variety of ways using various types of boundary conditions [|l|, ^, ||, |6|, |^, |8| . These calcu- 
lations, in general, rely on empirical potentials to simulate various configurations and incorporate 
a relatively large number of atoms. In general three types of boundary conditions are used, fixed 
cylindrical boundary conditions |2|, lattice greens functions |^ § or periodic boundary 

conditions 0]. For these approaches a large number of atoms is generally taken in order to mini- 
mize any artificial effects the boundary conditions may impose. Generally, this number is too large 
to be suited for ab initio calculations, even on the fastest and largest super-computers. Applying 
corrections to these effects is therefore important in order to have a reliable ab initio calculation. 

The simplest approach is to use fixed cylindrical boundary conditions. However, care must be 
taken due to the mismatch of the boundary with the dislocation, particularly when the dislocation 
moves. This mismatch can be minimized by performing calculations in very large systems, effec- 
tively extrapolating to the limit of infinite cylinder size. The size of the system needed to extract 
accurate values can be greatly reduced by applying leading order corrections due to the mismatch 
in the boundary |2|. However, even with these corrections, accurate calculations can generally be 
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made only for a cylinder with a radius of ~30A or greater, which corresponds to -^700 atoms when 
using periodic boundary conditions along the dislocation line. Such fixed cylindrical boundary 
conditions are generally ill-suited for an ab initio calculation due to both this large number of 
atoms and the artificial effects of the free surface at the boundary on the electrons. An alternative 
approach to account for boundary effects, and thereby reduce the size of the system needed for 
accurate calculation, is to use lattice greens function techniques |§|, § § • This technique, although 
quite elegant, still faces the issue of free surface effects when applied to first principles quantum 
mechanical calculations. 

Periodic boundary conditions are in general the most natural and straight forward way to cal- 
culate the Peierls stress in a density functional theory calculation. In periodic boundary conditions 
one is forced to have a net zero Burgers vector per unit cell. This is generally accomplished by using 



either a dipole or quadrupole array of dislocations 11, [1^. Peierls stress calculations have been 
done in periodic boundary conditions using empirical potentials, however, to avoid dislocation- 
dislocation interactions, generally requires the use of unit cells |Q] which would be impracticable 
for ab initio calculations. Wang et a/.[^ have introduced a promising method for use in small, 
periodic cells, but this approach requires the decomposition of the total energy into individual 
atomic contributions, which while feasible for inter-atomic potential approaches, is not well-defined 
in quantum mechanical calculations. 

In this paper it is shown how an accurate Peierls stress can be obtained in very small periodic 
cells. It is shown how leading order effects, due to the array of closely packed dislocations, can be 
accurately accounted. It is further shown how an accurate Peierls stress can be obtained by simply 
applying a pure shear to the system and minimizing the energy. This is a tremendous advantage in 
ab initio calculations as it minimizes the number of atomic configurations which need be explored. 
As an application of these ideas, we focus specifically on the calculation of the Peierls stress in bcc 
tantalum for a [lll]-screw dislocation in which the maximum resolved shear stress is along a {110}. 
This is among the most relevant geometries in understanding the plastic response of this material. 

The paper is organized as follows. Section 2 presents the results of Peierls stress calculations 
in very large cylindrical cells, as a reference against which we shall compare our new approach. 
Section 3 presents the results of calculations using periodic boundary conditions and underscores 
the problems that arise. Section 4 shows how to overcome these problems, and, finally. Section 5 
presents several techniques for accelerating the calculation. Finally, Section 6 concludes the paper. 



2 Reference Calculations 

All calculations presented throughout this manuscript employ a quantum-based, embedded atom 
force field (qEAMFF) developed in To provide the reference value of the Peierls stress for 

our calculations within periodic boundary conditions, we first perform calculations for isolated 
dislocations in a cylindrical geometry with fixed boundary conditions. 

To form the dislocation, we proceed as follows. First, we define coordinates for our calculations 
as follows: the x-axis lies along the [110] direction, the y lies along the [112], and the z lies along 
the [111]. A cylinder with radius R2 of Ta is then taken (Figure |l]) which has two regions, a fixed 
region and a relaxation region. In the relaxation region (r < Rl), where r is the radial distance 
from the center, the atoms are allowed to relax in response to their interatomic forces. The atoms 
in the fixed region {Rl < r < R2) are held fixed to their positions according to the solution of 
anisotropic elasticity theory. Finally, periodic boundary conditions (with period of one Burger's 
vector) are employed along the direction of the dislocation. Thus, these calculations consider only 
straight, infinite dislocations. A [111] -screw dislocation then is obtained by displacing all atoms in 
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the cyhnder according to the solution of anisotropic elasticity theory [14, ^] and then optimizing 
the positions of the atoms in the relaxation region. 

To calculate the Peierls stress, we then apply a strain to the system which ensures that the 
resulting stress has only one component, Oxz^ which generates a force on the dislocation line in the 
[112] direction [|l^. We then optimize the locations of the atoms in the relaxation region subject 
to this external strain. Figure shows the resulting dislocation structure in a large cyllinder 



(i?l ~ 1204) for a series of different strains using a differential displacement (DD) map [17|. In 
these maps, the circles represent columns of atoms viewed along the [111] direction. The arrows 
indicates the change in relative displacement that neighboring atomic columns make, relative to 
the bulk, due to the presence of the dislocation. The lengths of the arrows are normalized so 
that a displacement of 1/3 of a Burgers vector corresponds to an arrow of full length. The first 
panel (a) shows the ground-state structure of the dislocation at zero stress. The center of the 
dislocation is located at a diamond encased in triad of arrows. Going around this triad makes three 
displacements of one-third of a Burgers vector for a net displacement of one full Burgers vector 
relative to the bulk. The ground state structure of the core is seen to break the symmetry of the 
lattice by extending outward along three {110} planes. This "degenerate split core" is consistent 
with molecular dynamic results, using periodic boundary conditions, found in reference [|l| when 
using the same interatomic potential. 

As the strain is applied to the ground state structure (Figures ^-d), the dislocation feels force 
along the (llO)-plane. As the applied strain obtains a critical value, the dislocation center moves one 
lattice spacing along the (llO)-plane (Figure ^3). As the strain is increased further, the dislocation 
then glides along the (211) (Figures § c,d) in two steps, first along (101) and then spontaneously 
along (110). There are consequentially two Peierls barriers that the dislocation must overcome, 
the first along (110) and the second along (101). This leaves some ambiguity for the definition of 
the Peierls stress. Different authors have used different definitions [Q, ^, ^, ^]. For clarity we will 
consider the value of the axz stress for the first jump (Figure ^-b) as the first Peierls stress (PI) 
and the value of the Oxz stress for the second jump (Figure §b-c) as the second Peierls stress (-P2). 

To extract the limit of these critical stresses for an isolated dislocation in an infinite crystal, 
we have repeated the above calculations in cylinders of various sizes ranging from KV = 30^ to 
KV = I20A and extracted the Peierls stress as a function of radius. Figure ^ summarizes our results. 
In the limit, we find values of 0.74 GPa and 0.91 GPa for PI and P2, respectively. 

As a consistency check that we indeed expect finite-size effects to be small for our largest 
cylinders, one can employ the method of Shenoy and Phillips This methods estimates the 
unaccounted restoring stress that the boundary applies to a displaced dislocation to be 

Ksh , d 



where Ks = {Su/ {SuSii, — Si^))^^"^ in terms of the modified elastic compliances Sij [14], b is the 
Burgers vector, d is the distance the dislocation has moved from the center of the cylinder, R is the 
radius of the cylinder, and A is a dimensionless constant that can be calculated through elasticity 
theory or computationally. In our case, j4 ~ 2, Kg = 68GPa and b = 2.9A. Finally, for PI we 
have d ~ .lA and for P2 we have d ~ 2.7 A. For our largest cylinder (Rl = 120A), we calculate a 
restoring stress of 0.0004GPa and 0.012GPa, respectively. Both values are quite small, well within 
the uncertainties in Figure ^. 
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3 Periodic Boundary Conditions 



In periodic boundary conditions, the unit cell must contain a net zero Burgers vector. In practice 
this is generally accomplished through the use of either a dipolar or quadrupolar array of disloca- 



tions fl^, 11, 12 1 . Our calculations employ a quadrupole array, which has been shown to be the 



more appropriate choice for screw-dislocations |12|. Figure Q illustrates a cell containing 270 atoms. 

As our first attempt to calculate the Peierls stress in periodic boundary conditions, we take a 
quadrupole array with lattice vectors fixed at the values corresponding to the perfect bulk material, 
that is the appropriate lattice vectors prior the insertion of dislocations. Below, we refer to this 
choice as "unrelaxed" lattice vectors, as they generally do not correspond to the lattice vectors of 
an unstrained quadrupole array. 

To extract an estimate of the Peierls stress, we then follow a procedure analogous to that in 
Section ^, applying strain until the dislocations move. Here, rather than applying the strain to the 
fixed region (Rl < r < R2), we strain the unrelaxed lattice vectors. Then, we then compute the 
critical stress from the strain through the elastic constant matrix, 

a = C-e, (1) 

where o" is a column vector of the stresses, e is the applied strain vector, and C is the elastic 
constant matrix, assumed to equal that of the bulk material. The advantage of this approach is 
that it requires exploration of a minimal number of configurations and bulk lattice and elastic 
constants are relatively easy to obtain from first principles calculations. The disadvantage of this 
approach is that its underlying assumptions cast doubt of the accuracy of the results for unit cells 
of modest size. 

Figure ^ explores the convergence of this approach with increasing cell size for calculation of 
the first Peierls stress (-P1). The two smallest cells contain 90 atoms 2.9 A x 24A x 24A) and 
270 atoms (~ 2.9 A x 42^ x 42A) (or 45 and 135 atoms, if symmetry is exploited), and are the only 
cells convenient for detailed ab initio studies. The results in these cells, however, are extremely 
poor, with errors of 200% and 53%, respectively. Thus, great care must be taken when working 
with such small cells and a method is needed to correct for these finite-size effects. 



4 Corrections 

The preceding calculations make two major assumptions: (1) that the quadrupole array of dislo- 
cations does not change the elastic constant matrix from that of bulk material, and (2) that the 
use of unrelaxed lattice vectors is appropriate. To explore the impact of these assumptions, we 
repeat the above calculations while using both the relaxed lattice vectors and elastic constants of 
the quadrupolar array. 

Numerical calculations show that the symmetry of the elastic constant matrix of the quadrupolar 
array (C) is the same as that of the bulk, although the individual components may differ. In 
particular, the stress-strain relation takes the form 

a = C'-e (2) 
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In this equation, the lower 2x2 sub-block couples the cTxz and a^y stresses only to the corresponding 
strains, while the upper 4x4 block couples the crxx,o'yy,azz and ayz stresses only to their corre- 
sponding strains. Note, therefore, that for the present Peierls stress calculations only the lower 
2x2 block is relevant. 

Figure |^ shows the results (diamonds in the figure) of the extraction of the Peierls stress when 
using relaxed lattice vectors and the elastic constants for the quadrupole. The figure shows that the 
Peierls stress now converges much more quickly with cell size. (The figure does not include results 
for the smallest cell, which proved unstable to the relaxation of the lattice vectors.) For the smallest 
stable cell (135 atoms including symmetry), the error is reduced from 53% to only 18%, indicating 
the possibility of extracting reasonable results from cells of size suited to ab initio calculations. 
Although results for reasonably sized cells are accurate, this approach, however, is not necessarily 
well suited for ab initio calculations because relaxation of the lattice vectors and the calculation of 
the elastic constant matrix requires the exploration of many new atomic configurations. 



5 Accurate and Efficient Peierls Stress Calculation 

Having found an accurate approach, we now explore how to minimize that calculations associated 
with the above corrections. We begin by considering the necessity of relaxing the lattice vectors 
and then consider computation of the relevant components of the elastic constant matrix C. 

5.1 Benefit of fixed lattice vectors 



Working with unrelaxed lattice vectors generates spurious strains in the unit cell. The question 
is whether these strains create spurious stresses which confound the extraction of the Peierls stress. 
The elastic force which any such stresses would generate on the dislocation take the form p6[ |: 

FL = h-axr]. (4) 

Here, is the force per unit length on the dislocation, b is the Burgers vector, a is the stress tensor 
written as a matrix, and rj is the sense vector, the direction along which the dislocation runs. For 
a [lll]-screw dislocation b and tj both lie along the z— axis, and, therefore, only two components 
of stress can generate a force on the dislocation, axz and ayz, and confound the extraction of the 
Peierls stress. We now consider whether the strains associated with the unrelaxed lattice vectors 
can generate such stresses. We distinguish two types of strain, dilation and shear. 

The presence of the quadrupolar array tends to dilate the unit cell in the (111) plane. Due to 
symmetry, this dilation tends to be uniform, as we have confirmed by direct numerical calculation 
on the unit cells. The above form of the elastic constant matrix C' prevents such dilation from 
generating any a^z component. The second confounding component ayz would also vanish were 
the 4x4 subblock of the elastic constant matrix to have precisely the same symmetry as that of 
the bulk. Numerically, we find that that this is almost the case, and that the dilation contribution 
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to ayz ranges from 0.02% (in our largest cell) to only 4% (in the 135 atom cell) of the stress 
experienced in the unrelaxed cell. The diagonal components of strain en do not significantly affect 
the extraction of the Peierls stress, and therefore need not be relaxed to extract meaningful Peierls 
stresses. 

Turning to shear strains, the equal spacing of the dislocation array leads to zero shear within 
linear elasticity theory. From symmetry, however, core-core interaction (non-elastic) can generate 
a force in the [110]-direction only. Figure 0, which results in one non-zero component of strain, eyz- 
From direct calculations it is shown that this is in fact the only non-zero component of shear strain. 
Relaxing this component of strain, therefore, creates a artificial material environment, different 
from what isolated dislocations would experience. To more quickly approach the limit of isolated 
dislocations, one should therefore not relax this component of the strain. This advantage of using 
of unrelaxed lattice vectors is well known and commonly exploited in the materials literature in the 
study of iwo-dimensional defects such as grain boundaries or surfaces p^, |l9[. We therefore expect 
to be able to extract accurate Peierls stresses without the need for relaxing either off-diagonal or 
diagonal components of strain. 

Figure ^ shows the results (circles in the figure) of extracting the Peierls stress when using 
unrelaxed lattice vectors, but while still computing the stresses with the appropriate elastic constant 
matrix C. The results converge even more quickly than those obtained by relaxing the lattice 
vectors, thus supporting our analysis. Therefore, extremely good results (already within 2% in the 
135 atom cell) can be obtained by not relaxing the lattice vectors. This not only produces results 
of far superior quality, it also reduces the computational effort. 



5.2 Extraction of elastic constants 



The results of Section 5J., while quite impressive, still require calculation of the elastic constant 
matrix C of the quadrupolar array. Significant savings can be gained with the simple realization 
that only the 2x2 sub-block of the elastic constant matrix enters calculation of PI and P2 through 
Eq. ^. Moreover, within linear elasticity theory, these components can be extracted equally well 
at the relaxed or unrelaxed lattice vectors, thereby again mitigating the need for relaxation of the 
lattice vectors. 



Not all of the components of the 2x2 sub-block. 



44 ~*-^16 
16 *-^55 



(5) 



need be computed independently. 

To see this, consider the application of a pure, shear strain exz to the system. As a result of 
the symmetry of C, this generates only two components of stress, axz and axy From equation]^, 
however, axy does not generate any forces on the dislocation and thus, within linear elasticity, does 
not effect our calculation of the Peierls stress. The remaining stress, whose critical value is the 
Peierls stress, can be computed from just one component of C'2, 

= Cm^xz- (6) 

Finally, we note that this pre-factor can be extracted without any additional calculation. Figure ^ 
shows the total energy of the quadrupole array plotted as a function of strain during the extraction 
of the critical Peierls stress. Prior to the first dislocation glide event at exz ~ 0.015, the energy 
increases quadratically according to 

AE = \c',^elz, (7) 
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which contains precisely the same pre-factor as in Eq. ^. 

Figure ^ shows the result (triangles in the figure) of the extraction of the Peierls stress from the 
application of a pure, shear strain e^z to the unrelaxed lattice vecotrs and extracting the relevant 
elastic constant from the curvature of the energy prior to the glide events. The results are of nearly 
the same quality as working with unrelaxed lattice vectors and using the full elastic constant matrix 

c. 

5.3 Second Peierls stress: P2 

Figure ^ shows our preliminary results for the the extraction of the second Peierls stress P2. These 
results are complicated by the fact that after the first transition PI, the quadrupole array is 
now distorted and no longer a perfect quadrupole, thereby generating non-negligible dislocation- 
dislocation forces and further modifying the elastic constant matrix of the cell. As a possible 
correction to this effect, we are presently exploring moving the distorted cores back to their original 
quadrupole locations before further increasing the stress. Nonetheless, we find that significant 
improvements can be made by using the elastic constant matrix of a perfect quadrupole array and 
ignoring the dislocation-dislocation interactions. Finally, results for pure shear calculations are 
quite encouraging. 

6 Conclusion 

This paper present an accurate and effective way to calculate the {110} Peierls stress in a [111] 
screw dislocation in a bcc material. The results show that accurate results can be obtained even for 
the smallest cells while using unrelaxed lattice vectors and extracting the elastic constants directly 
from the calculations. The method most importantly appears to make ab initio calculations of 
Peierls stresses viable in periodic boundary conditions for the first time. 
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Figure 1: Pcierls stress calculation within cylindrical boundary conditions: atoms whose distance 
from the center is less than Rl are relaxed under inter-atomic potential forces, while those in the 
region between Rl and R2 are held fixed to the anisotropic elasticity theory solution. 
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Figure 2: Differential displacement (DD) maps of a [111] screw dislocation at increasing stress 
computed within cylindrical boundary conditions: a) zero applied stress, b) initial jump along 
(lIO) at a stress > 0.74 Gpa (PI), c) second jump, along (211) at a stress > 0.91 GPa (P2), d) 
subsequent jump along (211), e) orientation of {110} and {112} planes. 
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Figure 3: Convergence of PI and P2 as a function of 1/Rl, for cylinders ranging from Rl = 
to Rl = UOA. 
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Figure 4: Dislocation displacement (DD) map of a quadrupole dislocation array within a 270 atom 
cell with periodic boundary conditions. The cell may be reduced to 135 atoms with appropriate 
choice of lattice vectors. 
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Figure 5: Convergence of PI calculated within periodic boundary conditions using unrelaxed lattice 
vectors and bulk elastic constants: Peierls stress calculated using Eq. |l| (squares), limit exacted 
from Figure m (solid line with width indicating numerical uncertainty). Uncertainties in the periodic 
calculations are smaller than the square symbols. 
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Figure 6: Convergence of PI calculated within periodic boundary conditions: using relaxed lat- 
tice vectors and quadrupolar elastic constants (diamonds), using unrelaxed lattice vectors and 
quadrupolar elastic constants (circles), using unrelaxed lattice vectors, pure shear and elastic con- 
stants extracted during the calculation (triangles), asymptotic result extracted from cylindrical 
boundary conditions (horizontal line). Error bars associated with quadrupolar array are generally 
smaller than the associated icons. 
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Figure 7: The first nearest neighbors cannot generate any elastic force on the center dislocation, 
since they are equally spaced. Prom the core asymmetry there can be a core-core force in only in 
the [110]. Note that these symmetry arguments works for higher order neighbors too. 
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Figure 8: Energy of a 42A x A2A x 2.9A cell (270 atoms) plotted as a function of applied pure shear 
strain e^^: direct calculations (squares), quadratic fit (curve). 
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Figure 9: Convergence of P2 calculated within periodic boundary conditions: using unrelaxed lattice 
vectors and bulk elastic constants (squares), using relaxed lattice vectors and quadrupolar elastic 
constants (diamonds), using unrelaxed lattice vectors and quadrupolar elastic constants (circles), 
using unrelaxed lattice vectors, pure shear and elastic constants extracted during the calculation 
(triangles), asymptotic result extracted from cylindrical boundary conditions (horizontal line). 
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